Source code for hysop.symbolic.directional

# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.


import warnings
import numpy as np
import sympy as sm

from hysop.tools.htypes import first_not_None, check_instance
from hysop.symbolic import space_symbols
from hysop.symbolic.field import div, grad
from hysop.tools.warning import HysopWarning
from hysop.tools.sympy_utils import get_derivative_variables


[docs] def collect_direction(expr, var): is_tensor = isinstance(expr, np.ndarray) if is_tensor: R = np.empty_like(expr) for idx in np.ndindex(*R.shape): R[idx] = collect_direction(expr[idx], var) return R else: args = () if isinstance(expr, sm.Symbol): return expr elif isinstance(expr, sm.Derivative): if set(get_derivative_variables(expr)) != {var}: return 0 else: return expr elif isinstance(expr, (sm.Number, int, float)): return expr elif isinstance(expr, sm.Expr): for e in expr.args: de = collect_direction(e, var) args += (de,) return expr.func(*args) else: msg = ( f"Unknown expression type {type(expr).__mro__}.\nExpression was: {expr}" ) raise TypeError(msg)
[docs] def split(F, fixed_residue, force_residue, coords): from hysop.symbolic.relational import Assignment check_instance(F, sm.Basic) if isinstance(F, Assignment): msg = f"Expression cannot be of type Assignment: {F}" raise TypeError(msg) res = {} coords = first_not_None(coords, space_symbols) for i, xi in enumerate(coords): res[i] = collect_direction(F, xi) try: if force_residue is not None: residue = force_residue else: residue = (F - sum(res.values())).simplify() except AttributeError: msg = "Failed to compute residue for expression {}." msg = msg.format(F) raise RuntimeError(msg) count = sum((r != 0) for r in res.values()) if count > 0: for i in res.keys(): res[i] += residue / count if fixed_residue is not None: for i, r in res.items(): if r != 0: res[i] += fixed_residue return res
[docs] def split_assignement(expr, fixed_residue, force_residue, coords): lhs, rhs = expr.args res = split(rhs, fixed_residue, force_residue, coords=coords) for i, ei in res.items(): if ei != 0: res[i] = expr.func(lhs, ei) return res
if __name__ == "__main__": from hysop import Field, Box from hysop.tools.sympy_utils import enable_pretty_printing, greak from hysop.parameters.scalar_parameter import ScalarParameter enable_pretty_printing() box = Box(length=(1,) * 3) frame = box.frame U = Field(domain=box, name="U", is_vector=True) W = Field(domain=box, name="W", is_vector=True) alpha = ScalarParameter(name=greak[0]).s u = U.s(*frame.coords) w = W.s(*frame.coords) expr = (alpha * grad(u, frame) + (1 - alpha) * grad(u, frame).T).dot(w) # expr = div(np.outer(u,w), frame) print(expr) for i, x in split(expr, frame.coords).items(): print(i, x.simplify())